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Based on Brownian Dynamics computer simulations in two dimensions we investigate aggregation scenarios of colloidal particles 
with directional interactions induced by multiple external fields. To this end we propose a model which allows continuous change 
in the particle interactions from point-dipole-like to patchy-like (with four patches). We show that, as a result of this change, the 
non-equilibrium aggregation occurring at low densities and temperatures transforms from conventional diffusion-limited cluster 
aggregation (DLCA) to slippery DLCA involving rotating bonds; this is accompanied by a pronounced change of the underlying 
lattice structure of the aggregates from square-like to hexagonal ordering. Increasing the temperature we find a transformation to 
a fiuid phase, consistent with results of a simple mean-field density functional theory. 


1 Introduction 

Recent progress in the synthesis and directional binding of 
nanometer to micrometer sized patchy and anisotropic parti¬ 
cles makes possible the assembly of colloidal structures with 
multiple directed bonds^Sl^. The directional bonding can also 
be achieved by permanently embedded or field-induced dipole 
and/or multipole moments allowing directional and selective 
particle bondingHH. Within this class, particles with field- 
induced dipolar interactionsEEH are especially interesting be¬ 
cause switching the fields on and off is equivalent to switch¬ 
ing the particle intera ctions on and off. This means that ag¬ 
gregation mechanismstEMI ’dialed in’. Furthermore, 

the orientation of inducti ve fields may be used to direct par¬ 
ticle aggregation®32nil. In consequence, such directed self- 
assembly processes may be exploited for the formation of new 
functional materials with specific and/or adjustable proper¬ 
ties. Hence, understanding the interplay between externally 
induced particle properties, external fields and thermodynamic 
conditions, e.g., temperature, is of fundamental interest in 
modem material science, but also from a statistical physics 
point of view. 

An important subset of the many classes of self-assembled 
structures are percolated colloidal networks, which are char¬ 
acterized by system-spanning cross-linked (patchy) parti¬ 
cle cluste rs that are realizable even at low volume frac- 
tionsE^EStm, Such network-like aggregates are considered 
to be the underlying micro-structures of gels and have 
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been intensively investigated in experiment and theor y un- 
der equilibrium as well as non-equilibrium conditionsl ^^ l ^^ i 
In the latter, qualitatively different aggregation mechanisms 
can be identifi ed, n amely diffusion limited cluster aggre¬ 
gation (DLCX^^ and reaction limited cluster aggregation 
(RLCA)I^. In the DLCA regime each particle collision leads 
to the formation of a rigid and essentially (on the timescale of 
the experiment) unbreakable bond with fixed spatial orienta¬ 
tion. In contrast, in the RLCA regime the probability to form 
a rigid bond at collision is small. Systems with DLCA un¬ 
dergo irreversible dynamics and form fractal aggregates with 
specific fractal di mens ions Df « 1.71 in continuous two- 
dimensional spaceEI^. Such colloidal systems are consid¬ 
ered to be ’chemical gels’ and can be realized by having par¬ 
ticle interactions that are much stronger than ksT, preventing 
particles from dissociating due to thermal fluctuations. This 
leads to a pronounced hin drance of structural reconfiguration 
of large particle aggregates*^^*^. However, at higher tempera¬ 
tures these systems become ’physical gels’ where single par¬ 
ticles and larger substructures start to connect and disconnect 
fre quentl y. This strongly affects (increases) the fractal dimen- 
sionEmil and finally allows the system to achieve its equilib¬ 
rium state. 

A recently introduced new type of DLCA, which accounts 
for local rearrangements via fiexible bo nds, is slippery diffu¬ 
sion limited cluster aggregation (sDLCA)l^^^. Slippery bonds 
allow particles to move or rotate around each other as long 
as they stay in contact, meaning that bonds are still unbreak¬ 
able but can change their orientation. This additional degree 
of f reedo m generates, at least in three-dimensional simula- 
tionsl^S^, aggregates of the same fractal dimension as clas¬ 
sical DLCA but with a larger coordination number. 
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DLCA processes have been studi ed ext ensively in systems 
with isotropically attractive particles^^^^ but also in systems 
with patchy particles bearing permanent and/or locally re¬ 
stricted interaction sites on their surfaces®®. In the latter, 
the s patial or ientations of int eraction si tes can either be free to 
rotate®®® or fixed in spaced®®®. When the orientations 
of interaction sites are fixed in space, the associated ’chemical 
gels’ undergo anisotropic diffusion limited aggregation which 
yields a fractal dimension of Df ^ 1.5®®, lower than for the 
isotropic case. T his si tuation occurs, e.g., due to the presence 
of external fields®® or in lattice models®®, where motion 
is naturally restricted to certain directions. 

In the present paper we are particulary interested in the 
aggregation of colloids with field-induced multipolar inter¬ 
actions. Examples are capped ( metal -coated dielectric parti¬ 
cles studied earlier by some of where time-dependent 

electric fields can induce quadrupolar-like interactions. Here 
we consider even more complex interactions caused by 
crossed (orthogonal) fields. We briefly mention two examples 
of possible experimental realizations of such systems. The 
first one is a quasi two-dimensional system of suspended col¬ 
loidal particles, each composed of super-paramagnetic iron- 
oxide aggregates embedded in a polymer m atrix , which has 
been investigated experimentally by one of us^^®. In this case 
crossed external electric and magnetic fields, oriented in plane 
but perpendicular to each other, can be used to induce inde¬ 
pendent electric and magnetic dipole moments in the colloids 
leading to a directed self-assembly process resulting in two- 
dimensional single-particle chain networks. A second possi¬ 
ble experimental and quasi two-dimensional system consists 
of suspended colloidal particles under the infiuence of two in¬ 
plane orthogonal AC electric fields with a phase shift of tt. 
The fields will polarize the particles’ ionic layer periodically 
but at different times due to their phase shift. By adjusting the 
field frequencies and phases to the relevant timescales govern¬ 
ing particle diffusion and the relaxational dynamics of the po¬ 
larized ionic layer, two decoupled orthogonal dipole moments 
in each particle can in principle be generated by this setup. In 
both cases, the crossed dipole moments might be characterized 
as point-like or having a finite distance between their consti¬ 
tutive charges (or microscopic dipole moment distributions in 
the magnetic case). 

Here, we investigate the structure formation in such systems 
in a conceptional fashion by means of two-dimensional Brow¬ 
nian dynamics (BD) simulations of a generic particle model. 
The idea is to mimick externally-induced dipole moments via 
two pairs of screened Coulomb potentials that are decoupled 
to account either for magnetic and electric interactions or for 
two temporarily present electric interactions. The two charges 
associated with each pair are shifted outward from the par¬ 
ticle center, one parallel to the corresponding field and the 
other one anti-parallel. A sketch of such a particle with its 



Fig. 1 (Color online) Distribution of externally induced fictitious 
’’charges” q inside a particle. Positions of charges are determined by 
the vectors 5^^ G [—Sox^Sex, —Sey,Sey]) pointing either parallel 
or anti-parallel to the corresponding fields. 

internal arrangement of interaction centers is shown in Fig.[2 
By changing the charge separation, we systematically investi¬ 
gate the (transient) structural ordering and aggregation behav¬ 
ior predicted by this model. 

Highlights of our results are the following: At very high in¬ 
teraction energies and large charge separations we find that the 
particles undergo anisotropic diffusion limited cluster aggre¬ 
gation with rectangular local particle arrangements. Lowering 
the charge separation shifts the model behavior to a slippery 
diffusion limited aggregation (sDLCA) regime accompanied 
by a sharp transition of the lattice structure from rectangu¬ 
lar to hexagonal. In the proximity of this transition we ob¬ 
serve long-lived or arrested frustrated structures consisting of 
strongly interconnected hexagonal and rectangular lattice do¬ 
mains connected with each other. We also show that, upon in¬ 
crease of the temperature, the systems enter a fluid state. The 
corresponding ’fiuidization’ temperature turns out to be very 
close to the spinodal temperatures obtained from a mean-field 
density functional theory. 

The rest of this paper is oranized as follows. In section 
we present our model and discuss relevant target quantities ob¬ 
tained in the BD simulations. Numerical results are described 
in section!^ where we discuss first a specific low-temperature, 
low-density, state and then turn to the role of temperature and 
density. Finally, our conclusions are summarized in section]^ 

2 Theoretical Model 

We consider a two-dimensional system of N soft spheres of 
equal diameter a. The soft sphere interactions are repulsive 
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and are modeled by a shifted and truncated (12,6) Lennard- 
Jones Potential 


1 


Ussirij) = 4e - (cr/rijf + 1/4) (1) 

which is cut off at = 2^/®cr. Here, Vij = |rj — rj| is the 
particle center-to-center distance and e sets the unit of energy. 

The crossed orthogonal external fields induce orthogonal 
dipole moments = fiem and = jj.ee which we term 
for simplicity as ’magnetic’ and ’electric’ dipoles (although 
the model is also appropriate for two electric moments). The 
coordinate frame is adjusted to coincide with the directions 
of these moments so that = Bx and Oe = By. In gen¬ 
eral these moments could have different absolute values but 
for simplicity they are assumed to be equal. The two types 
of dipole moments are also assumed to be independent from 
each other and interact only with dipole moments of the same 
type on other particles. Intuitively, one would model the in¬ 
teraction energy between dipoles of particles 1 and 2 by the 
point-dipole potential 


TTOC o(A^ 1 •*'12)(A^2 •ri2) 

^dipi^ 12) - —3 3 -5 

'12 '12 


( 2 ) 


where a indicates the dipole type as being either e or m. Due 
to the constraint fii \\ /j 2 it follows that 




( 3 ) 
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The resulting total dipolar interaction between two particles is 
the sum of the dipolar potentials stemming from the magnetic 
and electric dipoles, respectively. Using jj = \fjf\ and the 
relation (ri 2 • Be + ri 2 • = rf2 (which holds since 

and are orthogonal) we obtain 

f/,%(ri2) + f/,™(ra2) = -^. (4) 

'12 

The resulting interaction on the right side of Eq. 0 is an 
isotropic, purely attractive interaction that lacks any kind 
of directional character. Therefore, the potential defined in 
Eq. Q can not gene rate any rectangular structures as observed 
in experimentsE2E2. Moreover, having in mind that the field- 
induced spatial separation of charges in the ionic layer of a 
suspended particle is comparable to the particle diameter cr, a 
point dipole model seems even more unreasonable. 

We therefore introduce an alternative model. Each dipole 
moment (with a = e, m) is replaced by two opposite 
charges — which are shifted out of the particle cen¬ 

ter by a vector 5^^ = (— l)^^eQ,, with /c = 1,2. The vec¬ 
tor 5"^^ points either parallel {k = 2) or antiparallel (k = 1) 
along the corresponding point dipole moment Indepen¬ 
dent of their shift or type, we set all charges to the same abso¬ 
lute value q = | = 2.5(e/cr)“^/^. The choice of the value 



Fig. 2 (Color online) Normalized direction-dependent pair interac¬ 
tion U(Yij) [see Eq. ([^] between a particle in the center of the co¬ 
ordinate frame and a second particle (indicated as black circle in 
(a)) at various positions Yij for three different charge separations 
6 = 0.1,0.21,0.3(7 corresponding to (a), (b), (c). Sterically ex¬ 
cluded areas are indicated by white circles, (d) Interaction energy 
at distance rij = cr as function of 0, the angle measured in multiples 
of TT against the x-axis, for 5 = O.lcr (yellow), 6 = 0.21cr (purple) 
and 6 = 0.3cr (black). 


2.5 is essentially arbitrary, as we will later normalize the in¬ 
teraction energy to eliminate the dependence of its magnitude 
on the charge separation 6 [see Eq. © below]. Charges k and 
I on different particles i and j interact via a Yukawa potential 


= -q^ 


r^OikOil 

ij 


(5) 


with = \rj — Vi — S^^\. The inverse screening 

length is choosen to /^ = 4.0cr“^ and a radial cutoff Tc = 4.0cr 
is applied. A schematic representation of the model with its in¬ 
ternal arrangement of ’charges’ is shown in Eig.[2 Due to the 
Yukawa-like interaction our model lacks the long-range char¬ 
acter of true dipolar interactions. However, similar models 
with comparable interaction ranges have been used to describe 
dipolar colloids in the c ontext of discontinous molecular dy¬ 
namics simulations®^. We note that mimicking magnetic 
dipoles via spatially separated ’charges’ appears somewhat ar¬ 
tificial from a physical point of view. However, here we could 
think of a particle with a strongly inhomogenous distribution 
of magnetic moments. Moreover, the idea behind our ansatz is 
not to mimick a particular complex colloid, but rather to pro¬ 
vide a generic model for field-induced directional interactions. 

The arrangment of charges inside particles then results in a 
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pair-interaction UDip{^ij) given by 

UDip{vii) = ^ (r,,) + (r,,)]. (6) 

k,l=l 


U = 7 = 1.0 is a friction constant and Ci(^*) is a 

Gaussian noise vector which acts on particle i and fulfills the 
relations (Ci) = 0 and (Ci(^)Ci(^O) = 2^kBT6ij6{t — t')® 
We perform simulations for up to lO^r^. 


In principle, UDipi^^ij) is a function of q and 5. To facilitate 
the comparison between the interactions at different 5 (q is 
chosen to be constant), we normalize UDip{^ij) according to 

Upipi^ij) = Uoipi^ij) X u/UDipicr^a) ( 7 ) 

where the constant u = —2.804e is calculated from the unnor¬ 
malized energy UDip{cF^a) with model parameters 5 = 0.3cr 
and q = 2.5(e/cr)“^/^. This procedure ensures that the nor¬ 
malized energy between two particles at contact (vij = a) and 
direction = crea (pointing along one of the fields) has the 
constant value u for all S, that is 

UDIp{cF^a) = ( 8 ) 


The full pair interaction of our model is then given by 


U{Yij) = Uss{^ij) + UDIp{^ij)- 


( 9 ) 


The resulting potential is illustrated in Fig.[^a)-(c) for a parti¬ 
cle in the center of the coordinate frame and a second particle 
at various distances and angles (j) = arccos(ri 2 • 
with ’charge’ separations 8 — 0.1,0.21,0.3cr. The value 
8 = 0 . 21(7 is motivated by our simulation results presented in 


Sec. |3.1| Sterically-excluded areas are shown in white and en¬ 
ergy values are color coded in units of e. The weak anisotropy 
of the resulting particle interactions at small 8 (where one es¬ 
sentially adds two dipolar potentials, see Eq. 0) transforms 
to a patchy-like pattern by increasing 8. Energy minima be¬ 
come more and more locally restricted and particle interac¬ 
tions become directional in character. This is also seen in 
Fig.j^d) which gives the energy between two particles in con¬ 
tact as function of 0 for different 8. From Fig. |^d) we also 
see that, independent of the ’charge’ separation 8, the min¬ 
ima of the full interaction potential [see Eq. 0 ] occur for 
connection vectors = crOe and (i.e., pointing 

along the fields). Note that this already holds for the unnor¬ 
malized energy given in. Eq. 0 . Simulations are performed 
with N = 1800 to 3200 particles at a range of reduced num¬ 
ber densities p* = pa^ and temperatures T* = /c^T/e, in a 
square-shaped simulation cell with periodic boundary condi¬ 
tions. The equations of motion 


N 




( 10 ) 


1=1 


are solved via the Euler scheme with an integration stepwidth 
At = 10 “^r 5 , where is the Brownian timescale defined by 


2.1 Target quantities 

To characterize the structure of the systems we consider sev¬ 
eral quantities. The first one is the mean coordination number 

1 ^ 
i=l 

where Zi is the number of neighbors of particle i and the sum 
is over all particles. In the following, two particles are consid¬ 
ered to be nearest neighbors if their center-to-center distance 
is smaller than = 1.15(7. 

To identify local particle arrangements, the orientational 
bond order parameter is of special importance. Eor particle 
k it is given by 

I (12) 

1=1 

with Zk being the number of neighbors and 0^^ = arccos(r/e/ • 
i*/ca/ being the angle between the bond of particle k 
and its neighbor I measured against a randomly chosen bond 
of particle k to one of its neighbouring particles A. Hence, 
= 0 for 2 (/c < 2. The integer value n determines the type 
of order which is detected by this parameter. We concentrate 
on 04 and 06 to identify square (rectangular) and hexagonal 
lattice types. Its ensemble average is calculated via 

1 ^ 

= (13) 

i=l 

The reversibility of ’bond’ formation and slipperyness of ex¬ 
isting bonds can be characterized by the bond and the bond- 
angle auto-correlation functions and Ca{t). To evaluate 
we assign a variable bij{t) to each pair of particles at 
each time step which is 1 if the particles i and j are nearest 
neighbors or zero otherwise. The bond auto-correlation func¬ 
tion is then defined as 

Cb{t) = {bij{to)bij{t)), (14) 

where the brackets indicate an average over all pairs that are 
bonded at time to. The bond-angle auto-correlation function 
Ca{t) is defined similarly by defining the unit vector 

=rij{t)/rij{t), (15) 

such that 

Ca{t) = (1 - arccos (aij(t) • aij{to))/7r) (16) 
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Fig. 3 (Color online) Simulation snapshots at p* = 0.3 and T* = 0.05 for (a) 6 = O.lcr, (b) S = 0.21cr and (c) S = 0.3cr. Particles are 
colored according to their value of 


where we average again over all pairs. While gives the 
information on how stable bonds are over time, Ca tells how 
stable their direction is over time. Note that in contrast to the 
typical definition of correlation functions for stationary sys- 
temsEil, here the functions and Ca{t) are not independent 
of the time origin to- 

Finally, we consider the fractal dimension Df of parti¬ 
cle clusters, which is particularly important in the context of 
DLCA. Clusters are defined as a set of particles with common 
next neighbors. The size of a cluster is then quantified by its 
radius of gyration 

Nci 

i=l 

where Nd is the number of particles in the cluster, and f is the 
postion of its center-of-mass. By plotting In Rg against In Nd 
for different clusters, we extract the fractal dimension Df via 
the relationship Rg ^ (see 


3 Results 


Our large-scale Brownian dynamics simulations show that the 
system is very sensitive to changes in temperature T*, num¬ 
ber density p*, and charge separation d. In this large parameter 
space we find a variety of different states ranging from small 
fractal aggregates and single-chain structures at low tempera¬ 
tures to coarser, isolated or interconnected clusters at higher 


temperatures. In the following sections 3.1 - 3.3 we first dis¬ 
cuss the structure, the time correlation functions and the frac¬ 
tal dimensions at a low temperature and an intermediate den¬ 
sity, focussing on the impact of the model parameter 5. In 
section |3.4| and |3.5| we then turn to the impact of temperature 
and density. 


3.1 Effect of Charge Separation on Local Order 

At first we study the system at low temperature T* =0.05 and 
intermediate density p* = 0.3 for different charge separations 
5. In Fig. [^simulation snapshots for 5 = O.la, 0.21cr, 0.3cr at 
t = 300r5 [see. Eq. below] are shown, where is the 
Brownian timescale. The colorcode reflects the orientational 
bond order parameter of each particle i. All three cases are 
characterized by clusters with irregular shapes. However, lo¬ 
cal particle arrangements differ strongly. While for 6 = 0.1a 
the particles aggregate in a hexagonal fashion, at J = 0.3cr 
they aggregate into rectangular structures. At the intermedi¬ 
ate charge separation 6 = 0.2la, hexagonal order dominates 
the system; however, some clusters also reveal subsets of par¬ 
ticles in rectangular arrangements. A more quantitative de¬ 
scription is given by the orientational bond order parameters 
4>4(e) shown in Fig. j^a) as functions of 6. By increasing 6, 
one observes a sharp transition at ^ 0.21a from hexagonal 

towards rectangular (square) order. 

The very presence of such a sharp transition can be ex¬ 
plained via energy arguments based on the (5-dependent pair 
potential plotted in Figs, [^a)-(d). To this end, we calculate 
the energy U{Yij) of a particle i with six 

neighbors j, which are located in a hexagonal arrangement at 
’contact’ distance a around i. Note that not all hexagonal con¬ 
figurations do have the same contact energy. This is due to the 
anisotropy of interactions, see Fig. [^d). Therefore we con¬ 
sider a hexagonal configuration in which the contact energy is 
as low as possible (this configuration was found numerically). 
The dependence of this lowest contact energy on the 

charge separation parameter is plotted in Fig. Also shown 
is the corresponding energy U-^{6) = ^= 4 x 

of a particle with four neighbors j located at distance a in a 
rectangular arrangement, i.e., in the energy minima around i 
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Fig. 4 (Color online) Results for simulations with N = 1800 at 
temperature T* = 0.05 and density p* = 0.3. (a) Orientational 
bond order parameters $4 for square (black) and $6 (yellow) for 
hexagonal particle arrangements, (b) Mean coordination number ^ 
as function of charge separation S at times t = 100, 200, SOOr?,. 


(the quantity u was defined below Eq. (|7])). Note that the en¬ 
ergy U-^{6) does not depend on 6 according to Eq. As 
shown in Eig. the two curves intersect at a ’’critical” value 
of (5 = 0.24cr. Thus, the simple energy arguments already 
suggest a transition between states with local hexagonal and 
square order, even though the predicted critical value is some¬ 
what larger than the value of = 0 . 2 la seen in the actual 
simulations at finite temperature and density [see Eig.j^a)]. 


Eurther information is gained from the behavior of the mean 
coordination number as a function of 6 plotted in Eig.^b) for 
three different times t = lOOr^, 200r5 and SOOr^. At all times 
considered, z undergoes a steep decrease aX 6 ^ 0 . 21 a from 
a nearly constant value, z^ex ~ 4.5, to a value Zsq ~ 3.5. 
This behavior reflects, on the one hand, again the presence 
of a sharp transition; on the other hand, the actual values of 
^hexi^sq) reveal the ”non-ideal” character of the aggregates in 
terms of coordination numbers. Eor example, for 6 > 0.21a 
we find that z and ^>4 decrease with S, while 4>6 increases. 
However, this does not indicate a decline of the rectangu¬ 
lar order; it rather results from an increasing amount of par¬ 
ticles residing in chains oriented either in x- or y-direction. 
The coordination number Zi of a particle i in such a chain is 
< 2, leading to a mean coordination number z < 4. Eurther- 
more, the parameters and [see Eq. (12)] become unity 



0.10 0.15 0.20 0.25 0.30 0.35 
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Fig. 5 (Color online) Minimum energy of a particle with six neigh¬ 
bors in hexagonal arrangement as function of S (black) and energy 
for a particle in rectangular arrangement with 4 neighbors (red). 


for a particle forming exactly two bonds under an angle of 
TT (straight chain). This does not affect 4 > 4 , which is already 
large at 6 > 0.21a, but significantly increases Finally, 
the counter-intuitive decrease of 4>4 with 6 results from the 
increasing amount of particles with only one neighbor (e.g., 
ends of chains appearing white in FigJ^c)). These particles 
yield no contribution to 4>4 [see Eq. (p^|)]. 

The ”non-ideal” values of Zhex and Zsq also explain why 
our energy argument for the location of the hexagonal-to- 
square transition, which was based on ideal arrangements with 
six and four neighbors, respectively, does yield the transition 
value S = 0.24a rather than 6 = 0 . 21 a obtained from sim¬ 
ulation. We can now reformulate the argument by using the 
actual mean coordination numbers extracted from our simu¬ 
lations, Zsq = 3.5 (instead of 4) and Zhex = 4.5 (instead 
of 6 ). Following the calculations for the ideal arrangements 
described before, the energy of the square-like arrangement is 
= 3.b XU. For the hexagonal arrangement, we use the av¬ 
erage minimum energy with either Zi = 4 or Zi = 5 neighbors, 
yielding U^{zhex^S) = {U^^^{4,6) There- 

suiting critical value of the charge separation is 6 « 0 . 21 a, 
which coincides nicely with the transition value observed in 
our simulations. 


3.2 Transient character of aggregates 

Although the local structures characterized by 2 and 4 > 4 (e) per¬ 
sist, in general, over the simulation times considered, we are 
still facing a transient (out-of-equilibrium) structure forma¬ 
tion as seen, e.g., from the slight increase of z with time in 
Fig.j^b). This raises a question about the typical ’’lifetime” of 
the aggregates. 

To this end we now consider dynamical properties, namely 
the bond and bond-angle auto-correlation functions, C})(t) and 
Ca{t). It is not reasonable to extract decay rates from these 
functions (as it is usually done) because in transient states, de¬ 
cay rates are, strictly speaking, functions of time themselves. 
Still, it is interesting to see whether the temporal correlation 
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Fig. 6 (Color online) Time correlation functions obtained from sim¬ 
ulations with N = 1800 at temperature T* = 0.05 and density 
p* = 0.3. (a) [(b)] Time evolution of the bond [angle] autocor¬ 
relation function Cb{t) [Ca{t)] for three different charge separations 
S = O.lcr, 0.21cr and 0.3cr colored in yellow, purple and black re¬ 
spectively. 


of bonds (bond angles) for different S allows us to distinguish 
between qualitatively different aggregation regimes. 

Numerical results for and Ca{t) are plotted in 

Figs.j^a) and (b), respectively, where we consider a large time 
range up to t lO^r^. The time axis starts at the finite time 
when all the systems have formed stable aggregates. The data 
in Figs, [^a) and (b) pertain to three representative values of 
the charge separation parameter related to the hexagonal struc¬ 
tures ((5 = O.lcr), rectangular structures ((5 = 0.3cr), and to the 
transition region ((5 = 0.21 cr). In the square regime ((5 = 0.3cr) 
the decay of both and Ca{t) is almost identical and very 
slow. From this we conclude that the square regime is charac¬ 
terized by almost unbreakable bonds with fixed orientations. 
This is different in the hexagonal regime (S = O.lcr) where 
remains nearly constant even after long times (meaning 
that bond-breaking is very unlikely), while Ca{t) decays much 
faster. Thus, the direction of bonds are less restricted. We 
interprete this behavior as evidence that two particles, though 
being bonded, are still able to rotate around each other to some 
extent. This is a characteristic feature of slippery bonds. Fi¬ 
nally, in the transition regime ((5 = 0.21cr) both functions Cb{t) 
and Ca{t) decay significantly faster than in the other cases, 
with the decay of the bond-angle correlation function being 
even more pronounced. In that sense we may consider the 
bonds in the transition region also as slippery (although less 
long-lived than in the other cases). 

We conclude that the different structural regimes identified 
in the preceding section are indeed characterized by different 
relaxational dynamics. Moreover, all of the observed aggre¬ 
gates have lifetimes of at least several hundered r^. Such long- 
lived bonds are indicative of diffusion limited cluster-cluster 
aggregation. In the next section we therefore consider the frac¬ 
tal dimension. 





Fig. 7 (Color online) Fractal dimension D/ as a function of charge 
separation at p* = 0.3 and T* = 0.05. At (5 = 0.21cr we find a 
bimodal distribution of fractal dimension with peaks at = 1.48 
(solid line) and 1.6 (dashed line). 


3.3 Diffusion limited aggregation 

In Fig. the fractal dimension Df\^ shown as a function of 
5 at time t = 250r5, density p* = 0.3 and temperature T* = 
0.05. We find that Df increases slightly with 6 but remains 
in a range between 1.4 and 1.5, except at ^ = 0.2la. There, 
the fractal dimension exhibits a bimodal distribution, taking 
values between Df ^ 1.48 and Df ^ 1.6 (dashed line in 
Fig.Q. 

Despite these variations, the values of Dp found here are 
significantly smaller than the fractal dimension Df = 1.71 ob¬ 
served in earlier studies of DL CA in two-dimensional contin¬ 
uous (off-lattice) systems^^^I^. Except for the case (5 = 0.21a, 
the values in Fig. are comparable with previous findings for 
DLCA in two-dimensional lattice systems and systems with 
spatial or interaction anisotropiesE^^E^. The present system is 
indeed anisotropic in the sense that the external fields impose 
preferences on the directions of particle bonds and therefore 
also on the orientations of aggregates. This effect is most 
pronounced in the rectangular regime ((5 = 0.3a). There¬ 
fore, it is plausible that our system undergoes a special case 
of anisotropic DLCA, in (quantitative) accord ance with ex¬ 
perimental results!^ and theoretical predictionsl^^I^^I^^I^. We 
should note that, due to our simulation method, the cluster 
sizes (typically involving 10^ — 10^ particles) are relatively 
small compared to th e partic le numbers considered in the lit¬ 
erature (10^ particles/^^J^^J^. A more detailed analysis of the 
impact of anisotropic interactions on the fractal structure is 
beyond the scope of this study. 

We also relate our findings to the newer concept of slip¬ 
pery DLCaI^^I^, where the bonds are essentially unbreakable 
but able to rotate. Indeed, as discussed in section [3^ bonds 
are slippery in nature for small d in the hexagona l reg ime. 
For three-dimensional systems it has been reportedl^^^H that 
the fractal dimension Df remains the same for slippery and 
classical DLCA, while the mean coordination number z dif- 
fers. Specifically, ^ is significantly higher for sDLCaE^SI^ 
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Fig. 8 (Color online) Temperature dependence of the system properties at density p* = 0.3 for charge separations S = 0.1,0.21,0.3cr colored 
in yellow, purple and black, respectively, (a) Fractal dimension Df evaluated ait ^ 250t 6, (b) Mean coordination number, (c) Orientational 
order parameters T >4 and T>6. 


The same observation emerges when we consider our values 
of z plotted in Fig. |^b), from which one sees a pronounced 
decrease of ^ upon entering the square (DLCA) regime. How¬ 
ever, in contrast to earlier studies we find to slightly in¬ 
crease with 6, especially in the hexagonal regime. We in¬ 
terpret this behavior as a consequence of the fact that bind¬ 
ing energies in the hexagonal regime decrease with increasing 
values of S, while they remain constant in the square regime 
(see Fig. |^. The corresponding stability of bonds should be 
correlated to the binding energies which explains the slightly 
increasing values of I)/ in the hexagonal regime. 

Finally, in the transition region ((5 = 0.2la) we found a bi- 
modal distribution of the fractal dimensions D f with maxima 
Sit Df « 1.48 and Df ^ 1.6. This second maximum corre¬ 
sponds to only « 25% of the considered cases (twelve inde¬ 
pendent simulation runs). The first maximum Df ^ 1.48 
therefore clearly dominates and fits nicely to the functional 
dependence of Df on S (see Fig. [^. We assume that the 
less frequent peak results from a switching of the local struc¬ 
tures between hexagonal and rectangular arrangements, which 
is accompanied by a significantly larger bond-breaking prob¬ 
ability (see Fig. |^c)). Again this allows compactification of 
aggregates and increases the fractal dimension in the transition 
regime. 

3.4 Beyond DLCA - Higher Temperatures 

Diffusion limited aggregation is restricted to systems with at¬ 
tractive particle interactions much stronger than ksT. By in¬ 
creasing the temperature sufficiently, thermal fiuctuations be¬ 
come able to break bonds which results in a faster decay of the 
the bond auto-correlation functions and a compactification of 
aggregates. Indeed, for square lattice models it was foun d that 
Dy is a monotonically increasing function of temperatureE3^. 

In Fig.j^a) the fractal dimension Df of the present model is 
plotted as a function of temperature T* for charge separations 
5 = 0.1a, 0.21a and 0.3a. 


We first concentrate on the case d = 0.3a, corresponding 
to the square regime at low T*. In the range of very low tem¬ 
peratures T* < 0.25, the fractal dimension is small and stays 
essentially constant. Increasing T* towards slightly larger val¬ 
ues then leads to an increase of Df, refiecting the (expected) 
compactification. This increase of Dy is accompanied by an 
increase of the mean coordination number ^ [see Fig. [^b)] 
within the temperature range considered, indicating the grow¬ 
ing number of bonds due to local and global structural re¬ 
configurations. The corresponding changes in the stability of 
the bonds are illustrated in Fig. where we have plotted the 
time evolution of C5(t) for several temperatures (at ^ = 0.3a). 
Clearly, the decay of becomes faster for higher tempera¬ 
tures. This is the reason why structural reconfigurations and, 
in consequence, compactification of aggregates becomes pos¬ 
sible. 

These trends persist until ^ 0.375, beyond which the 
system at S = 0.3a starts to behave in a qualitatively differ¬ 
ent way. The mean coodination number ^ displays a maxi¬ 
mum and subsequently a rapid decay. We also find that the 
fractal dimension has not yet reached its maximum value at 
this maximum occurs at the slightly larger temperature 
T* 0.42 (see Fig.[^a)). This ’delay’ of Dy can be under- 



Fig. 9 (Color online) Bond auto correlation function Cb{t) for dif¬ 
ferent temperatures T* at charge separation S = 0.3a and density 
p* = 0.3. 
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Fig. 10 (Color online) Simulation snapshots at p* = 0.3 with S = 0.21cr at (a) T* = 0.15, (b) T* = 0.3, and (c) T* = 0.45. Particles are 
colored according to their value (/)f. 


stood from the fact that, upon the entrance of bond-breaking, 
filigree parts of the aggregates are more likely affected than 
more compact ones. Hence, the fraction of ’compact’ small 
aggregates still grows. Even more important, the function 
^4(T*) in Fig. I^c) displays a pronounced decay of rectan¬ 
gular order for T* > Tf^gg. From the sum of these indications 
we conclude that, at temperatures higher than sq ^ 0.375, 
the system transforms into a (stable or metastable)phase. 
In this fiuid phase, the overall structure starts to become ho¬ 
mogeneous and isotropic, while the local structures involve 
only a small number of bonds with short bond-life times. 

For the system at = O.lcr (hexagonal structure at low T*), 
an estimate of the ’’fiuidization” temperature based on 

the behavior of order parameters, coordination number and 
fractal dimension is more speculative. Nevertheless, the data 
suggest that Tj5^. This is indicated, first, by the 
fact that decays only very slowly with temperature 

until T* 0.6 (see Fig.j^c)). Second, the mean coordination 
number shows only a weak maximum (and no fast decay after¬ 
wards) compared to the case 6 = 0.3cr. Third, the fractal di¬ 
mension keeps increasing with T* for all considered tempera¬ 
tures T* <0.6. Therefore we conclude that >0.6. We 

understand this higher fiuidization temperature at 6 = O.lcr 
from the fact that binding energies in hexagonal structures are 
larger; therefore, higher coupling energies must be overcome. 

To further justify these interpretations, particularly the 
emergence of fiuid phases, we performed a stability analy¬ 
sis of the homogenous isotropic high temperature state based 
on mean-field density functional theory (DFT). Specifically, 
we consider the isothermal compressibility xt- Positive 
values of xt imply that the homogeneous (fiuid) phase is 
stable, whereas negative values indicate that this phase is 
unstable. Specifically, the instability arises against long- 
wavelength density fluctuations, i.e. condensation. According 


to Kirkwood-Buff theoryl^ one has 


Xt^ oc 1 — pc{k = 0), (18) 

where c(0) is the Fourier transform of the direct correlation 
function (DCF) c(r 12) in the limit of long-wavelengths (k 
0). We approximate the DCF for distances Vij > a according 
to a mean field (MF) approximation, that is 

Cmf(i* 12) = —{kBT)~^U{Ti 2 ), ri2 > cr, (19) 

and use the Percus-Yevick DCF chs{'^ 12 ) of a pure hard- 
sphere fiuidl^for |ri2| < cr. The full DCF is then given by 

c(ri2) = CHs{ri2) + cmf(i*12)- (20) 

In Fig. m we present numerical results for the expression 
1 — pc{0) at p* = 0.3 as function of temperature. At low T*, 
all systems are characterized by negative values of 1 — pc(0). 
This indicates that the homogeneous isotropic phase is unsta¬ 
ble, consistent with the results of our simulations. Upon in¬ 
creasing T* the mean-field compressibility xt then becomes 



0.15 0.30 0.45 0.60 

* 


Fig. 11 (Color online) Numerical solutions to Eq 


18 


as function of 


T* for density p* = 0.3 and charge separations d — 0.1,0.21,0.3cr 
colored in yellow, purple and black, respectively. 
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indeed positive for all charge separations considered. Specif¬ 
ically, for 6 = 0.3cr the change of sign (related to a ’’spinodal 
point”) occurs at = 0.325 and for J = O.lcr at the much 
higher temperature = 0.6. These values are in surpris¬ 

ingly good agreement with our estimates for the ’’fluidization” 
temperatures based on the order parameter plots. 

The case S = 0.2la is again different. Here we And 
[see Fig.j^b)] that, starting from low temperatures inside the 
DLCA regime, the mean coordination number monotonically 
decreases. However, this does not indicate ’’fluidization” but 
rather a gradual transition from a state with dominant hexago¬ 
nal order towards a mixed state comprised of coexisting clus¬ 
ters with local hexagonal and square-like order. Indeed, [see 
Fig.j^c)], the orientational order parameters ^4 and ^6 reveal 
that the fraction of particles bound in square clusters increases 
with T* and Anally overtakes the fraction of particles involved 
in hexagonal clusters at T* ^ 0.35. Corresponding snapshots 


of simulation results are shown in Fig. 10 At all temperatures 
considered one observes separated clusters. With increasing 
temperature their shape becomes more regular, while the lo¬ 
cal rectangular order becomes more pronounced. Finally, at 
T* = 0.45 the fractal dimension Df and the square order 
parameter 4>4 reach their maximum values, suggesting a ’’flu¬ 
idization” similar to the behavior observed at other values of 5. 
Interestingly, our stability analysis [see Eq. ( p^ ] indicates an 
instability at the same temperature = 0.45. With this sur¬ 
prisingly accurate agreement between theory and simulation, 
we conclude that in the transition regime (^ = 0.21a), increas¬ 
ing thermal fluctuations first push the system from a domi¬ 
nantly hexagonal state into a rectangular one, which then en¬ 
ters a metastable fluid phase after passing the ’’spinodal point”. 


3.5 Spotlight on higher densities 

In this section we revisit the system behavior at the low tem¬ 
perature T* = 0.05, but consider different densities in the 
range p* < 0.7. Whereas low-density systems at T* = 0.05 
display DLCA as discussed in sections A-C, this aggrega¬ 
tion mechanism is expected to disappear at higher densities: 
here, the particles are just unable to diffuse sufficiently freely. 
Rather, the particles will very frequently collide and then im¬ 
mediately form rigid bonds. A typical structure at the high¬ 
est density considered, p* = 0.7, and separation parameter 


S = 0.21a is shown in Fig. 12 Clearly, the system is per¬ 
colated, that is, the particles form a single, system-spanning 
cluster. Interestingly, this cluster is composed of extended re¬ 
gions characterized by either square-like order or hexagonal 
order. We note that, at 6 = 0.21a, simultaneous appearance 
of clusters with both types of order also occurs at low den¬ 
sities and higher temperatures (see section [3^ . However, at 
the high density considered here the regions of each type are 
larger and the particle arrangements are much more regular 



Fig. 12 (Color online) System at T* = 0.05, density p* = 0.7 
and S = 0.21a. The colorcode gives the orientational bond-order 
parameter of each particle i. 


(i.e., there are less defects). 

To better understand the impact of the density on the cluster 
structures we plot in Fig. [^a) the orientational bond order 
parameters 4>4 and 4>6 as functions of p* for ^ = 0.21a (at 
6 = 0.1a and S = 0.3a the order parameters are essentially 
independent of the density). From Fig. [^a) it is seen that 
the amount of rectangular (hexagonal) order sharply increases 
(decreases) at a density of p* ^ 0.45. This is a surprising 
result as one would expect that, upon compressing the system, 
close-packed, hexagonal structures rather become more likely. 
However, at the low temperature considered here, structural 
reorganization is strongly hindered. 

We also note that all of the systems investigated at densi¬ 
ties p* > 0.45 turned out to be percolated (suggesting that the 
value p* = 0.45 is indeed related to the percolation transi¬ 
tion). It thus seems that the percolation tends to stabilize the 
initially formed square-lattice symmetry, as the subsequent re¬ 
organization is hindered by the lack of mobility. In effect, we 
are faced with quenched states that could not densify within 
the time domain studied. This interpretation is also consistent 
with the decrease of the mean coordination number once the 
system is percolated (p* > 0.45) as shown in Fig.p^b). 

4 Conclusions 

In this work we propose a new model for field-directed ag¬ 
gregation of colloidal particles with anisotropic interactions 
induced by external fields. Th e model was inspired by recent 
experimental workEEGIEEll on novel colloidal particles in 
which external fields can induce two, essentially decoupled, 
dipoles. 
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Fig. 13 (Color online) (a) Orientational bond order parameter $4 
and $6 as function of density p* for (5 = 0.21cr at T* = 0.05. (b) 
Mean coordination number z as function of p* at T* = 0.05 for 
different 6 = O.la, 0.21cr and 0.3cr colored in yellow, purple and 
black, respectively. 

The formation of particle networks with multiple percola¬ 
tion directions can find application in a range of new materi¬ 
als with anisotropic electrical and thermal conduction, mag¬ 
netic or electric polarizability or unusual rheological proper¬ 
ties. The aggregated clusters can be dispersed in liquid, while 
the percolated networks can be embedded in a polymer or gel 
medium® The key to the fabrication of such novel classes 
of materials containing particle clusters and networks is the 
control of the process parameters to obtain the desired inter¬ 
connectivity, density and structure. 

Against this background, the focus of our theoretical study 
was to understand the formation of transient, aggregated struc¬ 
tures appearing at low-temperatures. Performing large-scale 
BD simulations we have found that, depending on the distri¬ 
bution of the field-induced attractive ’’sites” in the particles, 
different aggregation mechanisms arise. These have been an¬ 
alyzed via appropriate structural order parameters, bond time- 
correlation functions as well as by the fractal dimension. Our 
BD results demonstrate that by varying the charge separation 
parameter, that is, the distribution of attractive sites, the sys¬ 
tems transform from DLCA (essentially rigid bonds) towards 
sDLCA (slippery bonds). Moreover, we show that the change 
of aggregation behavior is accompanied by significant changes 
of the local cluster structure. 

Indeed, the cluster structure can be easily manipulated by 
exploiting the interplay between temperature, density and 
model parameter d. This allows formation of unexpected 
structures e.g., pronounced rectangular packing instead of 
closed packed hexagonal structures by increasing density. 
This unusual behavior appears to be dictated by the inabil¬ 
ity of the originally formed lattices with square symmetry 
to re-arrange into more dense hexagonal lattices. It has po¬ 
tentially important consequences for colloidal assembly, as is 
points out the ability to use mutidirectional field-driven assem¬ 
bly for the making of lower-density, yet highly interconnected, 
phases. 


Future research should focus on a more detailed investiga¬ 
tion of the interplay between the aggregation mechanisms ob¬ 
served here (anisotropic and slippery DLCA), and the equi¬ 
librium phase behavior, particularly the location of a conden¬ 
sation transition and of percolation at higher densities. This 
includes investigation of the influence of entropy which we 
did not discuss but is expected to strongly influence the aggre¬ 
gation behavioi®. 

In conclusion, our study makes an important contribution to 
current research on structural and dynamical pheno mena ac¬ 
companying self-assembly of complex colloids^^^^^. In partic¬ 
ular, we have introduced a generic model describing colloids 
in multidirectional fields yielding tunable multipolar interac¬ 
tions. Our study thus provides a microscopic understanding of 
aggregation processes in such systems. 
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